#include <iostream>
#include<stdio.h>
#include<time.h>
#include<pthread.h>
using namespace std;
const int NUM_THREADS = 7;//工作线程数量
using namespace std;
const int N = 1000;
float m[N][N];
void m_reset() {
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < i; j++)
            m[i][j] = 0;
        m[i][i] = 1.0;
        for (int j = i + 1; j < N; j++)
            m[i][j] = rand() ;
    }
    for (int k = 0; k < N; k++)
        for (int i = k + 1; i < N; i++)
            for (int j = 0; j < N; j++)
                m[i][j] += m[k][j];
}
void serial() {
    for (int k = 0; k < N; k++)
    {
        for (int j = k + 1; j < N; j++)
            m[k][j] = m[k][j] / m[k][k];
        m[k][k] = 1.0;
        for (int i = k + 1; i < N; i++)
        {
            for (int j = k + 1; j < N; j++)
                m[i][j] = m[i][j] - m[i][k] * m[k][j];
            m[i][k] = 0;
        }
    }

}
typedef struct {
	int t_id; //线程 id
}threadParam_t; 
//barrier 定义
pthread_barrier_t barrier_Divsion;
pthread_barrier_t barrier_Elimination;
void* threadFunc(void* param) {
	threadParam_t *p = (threadParam_t*)param;
	int t_id = p -> t_id; 
	for (int k = 0; k < N; ++k) {
		// t_id 为 0 的线程做除法操作，其它工作线程先等待
		// 这里只采用了一个工作线程负责除法操作，同学们可以尝试采用多个工作线程完成除法操作
		if (t_id == 0)
			for (int j = k + 1; j < N; j++) {
				m[k][j] = m[k][j] / m[k][k];
			}
		m[k][k] = 1.0;


		//第一个同步点
		pthread_barrier_wait(&barrier_Divsion);
		//循环划分任务（同学们可以尝试多种任务划分方式）
		for (int i = k + 1 + t_id; i < N; i += NUM_THREADS) {
			//消去
			for (int j = k + 1; j < N; ++j) {
				m[i][j] = m[i][j] - m[i][k] * m[k][j];
			}
			m[i][k] = 0.0;
		}

		// 第二个同步点
		pthread_barrier_wait(&barrier_Elimination);
	}
	pthread_exit(NULL);
	return NULL;
}

int main() {
    struct timespec sts,ets;
    time_t dsec;
    long dnsec;
    
	m_reset();
	//初始化 barrier
	pthread_barrier_init(&barrier_Divsion, NULL, NUM_THREADS);
	pthread_barrier_init(&barrier_Elimination, NULL, NUM_THREADS);

    timespec_get(&sts,TIME_UTC);
    serial();
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0){
             dsec--;
    dnsec+=1000000000ll;
    }
    cout<<"LU分解单线程:";
    printf ("%lld.%09llds\n",dsec,dnsec);

    timespec_get(&sts,TIME_UTC);
   //创建线程
	pthread_t handles[NUM_THREADS];// 创 建 对 应 的 Handle 
	threadParam_t param[NUM_THREADS];// 创建对应的线程数据结构
	for (int t_id = 0; t_id < NUM_THREADS; t_id++){
		param[t_id].t_id = t_id; pthread_create(&handles[t_id], NULL, threadFunc, (void*)&param[t_id]);
	}

	for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
		pthread_join(handles[t_id], NULL);
	}

		//销毁所有的 barrier
	pthread_barrier_destroy(&barrier_Divsion);
	pthread_barrier_destroy(&barrier_Elimination);
    timespec_get(&ets,TIME_UTC);
    dsec=ets.tv_sec-sts.tv_sec;
    dnsec=ets.tv_nsec-sts.tv_nsec;
    if (dnsec<0){
             dsec--;
    dnsec+=1000000000ll;
    }
    cout<<"LU分解多线程:";
    printf ("%lld.%09llds\n",dsec,dnsec);

	
	

		return 0;
}